#!/usr/bin/env python3

import argparse
import xarray as xr
from multiprocessing import Pool
from itertools import zip_longest, repeat
import metpy.calc as calc
from metpy.units import units
import numpy as np
import pandas as pd
import pendulum
from scipy.interpolate import interp1d
from glob import glob
import os
import sys
sys.path.append(f'{os.path.dirname(os.path.realpath(__file__))}/../../utils')
from utils import cli, parse_time, parse_config, create_kdtree, find_bbox, create_bilinear_wgts, apply_bilinear, has_key

man_plev = np.array((1000, 925, 850, 700, 500, 400, 300, 250, 200, 150, 100, 70, 50, 30, 20, 10, 7, 5, 3, 2, 1)) * 100 * units.Pa

def sim_sound_ecwf(file_root, sid, lon, lat, plev, start_time, end_time):
	cli.notice(f'Simulate SOUND at {sid}@({lon},{lat}).')
	files = file_root + '/W_NAFP_C_ECMF_*'
	used_files = []
	# Assume the file name contains date time.
	times = []
	for file in sorted(glob(files)):
		tmp = os.path.basename(file)
		time = pendulum.datetime(year=int(tmp[14:18]), month=int(tmp[42:44]), day=int(tmp[44:46]), hour=int(tmp[46:48]))
		# f = xr.open_dataset(file, engine='cfgrib', backend_kwargs={
		# 	'filter_by_keys': {
		# 		'typeOfLevel': 'surface',
		# 		'edition': 1
		# 	},
		# 	'indexpath': ''
		# })
		# time = pendulum.parse(str(f.variables['time'].values))
		if start_time <= time <= end_time:
			times.append(time)
			used_files.append(file)
	if len(used_files) == 0:
		cli.error(f'Cannot find valid ECWF files!')
	f_sfc = xr.open_mfdataset(used_files, engine='cfgrib', backend_kwargs={
		'filter_by_keys': {
			'typeOfLevel': 'surface',
			'edition': 1
		},
		'indexpath': ''
	}, concat_dim='valid_time', combine='nested')
	f_plev = xr.open_mfdataset(used_files, engine='cfgrib', backend_kwargs={
		'filter_by_keys': {
			'typeOfLevel': 'isobaricInhPa',
			'edition': 1
		},
		'indexpath': ''
	}, concat_dim='valid_time', combine='nested')

	# Just use the pressure levels in the data
	# if type(plev) != list: plev = man_plev

	# Surface
	u10  = f_sfc.u10.interp(longitude=lon, latitude=lat, method='linear') * units.mps
	v10  = f_sfc.v10.interp(longitude=lon, latitude=lat, method='linear') * units.mps
	t2   = f_sfc.t2m.interp(longitude=lon, latitude=lat, method='linear') * units.K
	d2   = f_sfc.d2m.interp(longitude=lon, latitude=lat, method='linear') * units.K
	ps   = f_sfc.sp .interp(longitude=lon, latitude=lat, method='linear') * units.Pa
	slp  = f_sfc.msl.interp(longitude=lon, latitude=lat, method='linear') * units.Pa
	sh2  = calc.specific_humidity_from_dewpoint(ps, d2)
	mr2  = calc.mixing_ratio_from_specific_humidity(sh2)
	# tv   = calc.virtual_temperature_from_dewpoint(ps, t2, d2)
	tv   = calc.virtual_temperature(t2, mr2)
	gam  = 0.0065
	rd   = 287.05
	g    = 9.80665
	zs   = (float(tv) / gam * ((slp / ps)**(gam * rd / g) - 1)) * units.m
	ws10 = calc.wind_speed(u10, v10)
	wd10 = calc.wind_direction(u10, v10)
	rh2  = calc.relative_humidity_from_dewpoint(t2, d2) * 100 * units.percent
	sfc_data = {
		'sid' : [sid for _ in range(len(times))],
		'lon' : [lon for _ in range(len(times))],
		'lat' : [lat for _ in range(len(times))],
		'zs'  : [float(zs)],
		'time': times,
		'u10' : [float(u10)],
		'v10' : [float(v10)],
		't2'  : [float(t2)],
		'sh2' : [float(sh2)],
		'ps'  : [float(ps)],
		'slp' : [float(slp)],
		'ws10': [float(ws10)],
		'wd10': [float(wd10)],
		'rh2' : [float(rh2)],
		'd2'  : [float(d2)]
	}

	# Pressure levels
	plev = f_plev.isobaricInhPa.values * 100 * units.Pa
	h    = f_plev.gh.interp(longitude=lon, latitude=lat, method='linear').values / g * units.m
	u    = f_plev.u .interp(longitude=lon, latitude=lat, method='linear').values * units.mps
	v    = f_plev.v .interp(longitude=lon, latitude=lat, method='linear').values * units.mps
	t    = f_plev.t .interp(longitude=lon, latitude=lat, method='linear').values * units.K
	sh   = f_plev.q .interp(longitude=lon, latitude=lat, method='linear').values * units('kg/kg')
	rh   = calc.relative_humidity_from_specific_humidity(plev, t, sh).to(units.percent)
	d    = calc.dewpoint_from_specific_humidity(plev, t, sh).to(units.K)
	ws   = calc.wind_speed(u, v)
	wd   = calc.wind_direction(u, v)
	plev_data = []
	for i in range(len(times)):
		for j in range(len(plev)):
			if not np.isnan(u[i,j]):
				plev_data.append({
					'sid' : sid,
					'lon' : lon,
					'lat' : lat,
					'time': times[i],
					'p'   : plev.magnitude[j],
					'h'   : h .magnitude[i,j],
					'u'   : u .magnitude[i,j],
					'v'   : v .magnitude[i,j],
					't'   : t .magnitude[i,j],
					'sh'  : sh.magnitude[i,j],
					'd'   : d .magnitude[i,j],
					'rh'  : rh.magnitude[i,j],
					'ws'  : ws.magnitude[i,j],
					'wd'  : wd.magnitude[i,j]
				})
	return [pd.DataFrame(sfc_data), pd.DataFrame(plev_data)]

def sim_sound_wrf(file_root, sid, lon, lat, plev, start_time, end_time):
	cli.notice(f'Simulate SOUND at {sid}@({lon},{lat}).')
	files = file_root + '/wrfout*'
	used_files = []
	# Assume the file name contains date time.
	times = []
	for file in sorted(glob(files)):
		f = xr.open_dataset(file, engine='netcdf4')
		time = pendulum.from_format(str(f.Times[0].data)[2:-1], 'YYYY-MM-DD_HH:mm:SS')
		if start_time <= time <= end_time:
			times.append(time)
			used_files.append(file)
	if len(used_files) == 0: cli.error('Cannot find wrfout files within time range!')
	f = xr.open_mfdataset(used_files, concat_dim='Time', parallel=True, combine='nested', chunks={'south_north': 40, 'west_east':40, 'south_north_stag':40, 'west_east_stag':40})

	if type(plev) != list: plev = man_plev

	xlon = f.XLONG[0,:,:]
	xlat = f.XLAT [0,:,:]
	xlon_u = f.XLONG_U[0,:,:]
	xlat_u = f.XLAT_U [0,:,:]
	xlon_v = f.XLONG_V[0,:,:]
	xlat_v = f.XLAT_V [0,:,:]
	kdtree = create_kdtree(xlon, xlat)
	kdtree_u = create_kdtree(xlon_u, xlat_u)
	kdtree_v = create_kdtree(xlon_v, xlat_v)
	bbox = find_bbox(kdtree, lon, lat)
	bbox_u = find_bbox(kdtree_u, lon, lat)
	bbox_v = find_bbox(kdtree_v, lon, lat)
	wgts_h = create_bilinear_wgts(xlon, xlat, bbox, lon, lat)
	wgts_h_u = create_bilinear_wgts(xlon_u, xlat_u, bbox_u, lon, lat)
	wgts_h_v = create_bilinear_wgts(xlon_v, xlat_v, bbox_v, lon, lat)

	# Surface
	zs   = apply_bilinear(wgts_h, f.HGT , dims=('south_north','west_east')) * units.m
	u10  = apply_bilinear(wgts_h, f.U10 , dims=('south_north','west_east')) * units.mps
	v10  = apply_bilinear(wgts_h, f.V10 , dims=('south_north','west_east')) * units.mps
	t2   = apply_bilinear(wgts_h, f.T2  , dims=('south_north','west_east')) * units.K
	mr2  = apply_bilinear(wgts_h, f.Q2  , dims=('south_north','west_east'))
	ps   = apply_bilinear(wgts_h, f.PSFC, dims=('south_north','west_east')) * units.Pa
	slp  = ps * (1 + 0.0065 * units.K / units.m * zs / t2) ** 5.257
	ws10 = calc.wind_speed(u10, v10)
	wd10 = calc.wind_direction(u10, v10)
	rh2  = calc.relative_humidity_from_mixing_ratio(ps, t2, mr2).to(units.percent)
	d2   = calc.dewpoint_from_relative_humidity(t2, rh2).to(units.K)
	sh2  = calc.specific_humidity_from_mixing_ratio(mr2)
	sfc_data = {
		'sid' : [sid for _ in range(len(times))],
		'lon' : [lon for _ in range(len(times))],
		'lat' : [lat for _ in range(len(times))],
		'zs'  : list(zs.data),
		'time': times,
		'u10' : list(u10.data),
		'v10' : list(v10.data),
		't2'  : list(t2.data),
		'sh2' : list(sh2.data),
		'ps'  : list(ps.data),
		'slp' : list(slp.data),
		'ws10': list(ws10.data),
		'wd10': list(wd10.data),
		'rh2' : list(rh2.data),
		'd2'  : list(d2.data),
	}

	# Pressure levels
	# First, do horizontal interpolation.
	p      = apply_bilinear(wgts_h, f.P + f.PB, dims=('south_north','west_east'), debug=True) * units.Pa
	ph     = 0.5 * (f.PH [:,0:f.PH .shape[1]-1,:,:] + f.PH [:,1:f.PH .shape[1],:,:]) / 9.8
	phb    = 0.5 * (f.PHB[:,0:f.PHB.shape[1]-1,:,:] + f.PHB[:,1:f.PHB.shape[1],:,:]) / 9.8
	tmp_h  = apply_bilinear(wgts_h, ph + phb, dims=('south_north', 'west_east'))
	tmp_u  = apply_bilinear(wgts_h_u, f.U, dims=('south_north', 'west_east_stag'))
	tmp_v  = apply_bilinear(wgts_h_v, f.V, dims=('south_north_stag', 'west_east'))
	tmp_pt = apply_bilinear(wgts_h, f.T + 300, dims=('south_north', 'west_east'))
	tmp_mr = apply_bilinear(wgts_h, f.QVAPOR, dims=('south_north', 'west_east'))
	# Second, do vertical interpolation.
	h = []; u = []; v = []; pt = []; mr = []
	for i in range(p.shape[0]): # time dimension
		h .append(interp1d(p[i,:], tmp_h [i,:], kind='linear', bounds_error=False, copy=False)(plev))
		u .append(interp1d(p[i,:], tmp_u [i,:], kind='linear', bounds_error=False, copy=False)(plev))
		v .append(interp1d(p[i,:], tmp_v [i,:], kind='linear', bounds_error=False, copy=False)(plev))
		pt.append(interp1d(p[i,:], tmp_pt[i,:], kind='linear', bounds_error=False, copy=False)(plev))
		mr.append(interp1d(p[i,:], tmp_mr[i,:], kind='linear', bounds_error=False, copy=False)(plev))
	h  = np.array(h ) * units.m
	u  = np.array(u ) * units.mps
	v  = np.array(v ) * units.mps
	pt = np.array(pt) * units.K
	mr = np.array(mr) * units('kg/kg')
	# Third, calculate other quantities.
	ws = calc.wind_speed(u, v)
	wd = calc.wind_direction(u, v)
	t  = calc.temperature_from_potential_temperature(plev, pt)
	rh = calc.relative_humidity_from_mixing_ratio(plev, t, mr).to(units.percent)
	sh = calc.specific_humidity_from_mixing_ratio(mr)
	d  = calc.dewpoint_from_relative_humidity(t, rh).to(units.K)
	plev_data = []
	for i in range(len(times)):
		for j in range(len(plev)):
			if not np.isnan(u[i,j]):
				plev_data.append({
					'sid' : sid,
					'lon' : lon,
					'lat' : lat,
					'time': times[i],
					'p'   : plev[j].magnitude,
					'h'   : h [i,j].magnitude,
					'u'   : u [i,j].magnitude,
					'v'   : v [i,j].magnitude,
					't'   : t [i,j].magnitude,
					'sh'  : sh[i,j].magnitude,
					'd'   : d [i,j].magnitude,
					'rh'  : rh[i,j].magnitude,
					'ws'  : ws[i,j].magnitude,
					'wd'  : wd[i,j].magnitude,
				})
	return [pd.DataFrame(sfc_data), pd.DataFrame(plev_data)]

def add_sfc_obs_err(df, obs_err_u, obs_err_v, obs_err_t, obs_err_q, obs_err_p):
	df['u10' ] += np.random.normal(loc=0, scale=obs_err_u, size=len(df))
	df['v10' ] += np.random.normal(loc=0, scale=obs_err_v, size=len(df))
	df['t2'  ] += np.random.normal(loc=0, scale=obs_err_t, size=len(df))
	df['ps'  ] += np.random.normal(loc=0, scale=obs_err_p, size=len(df))
	for i in range(len(df)):
		sh0 = df.loc[i,'sh2']
		while True:
			sh = sh0 + np.random.normal(loc=0, scale=obs_err_q, size=1)[0]
			if sh < 0: continue
			rh = calc.relative_humidity_from_specific_humidity(df['ps'].iloc[i] * units.Pa, df['t2'].iloc[i] * units.K, sh)
			if rh > 1: continue
			break
		df.loc[i,'sh2'] = sh
		df.loc[i,'rh2'] = rh.magnitude * 100
	df['d2'  ] = calc.dewpoint_from_relative_humidity(list(df['t2']) * units.K, list(df['rh2']) * units.percent).to(units.K).magnitude
	df['slp' ] = df['ps'] * (1 + 0.0065 * df['zs'] / df['t2']) ** 5.257
	df['ws10'] = calc.wind_speed(list(df['u10']) * units.mps, list(df['v10']) * units.mps)
	df['wd10'] = calc.wind_direction(list(df['u10']) * units.mps, list(df['v10']) * units.mps)
	return df

def add_plev_obs_err(df, plev_obs_err, obs_err_u, obs_err_v, obs_err_t, obs_err_q):
	tmp = pd.DataFrame({
		'p': plev_obs_err,
		'u': obs_err_u,
		'v': obs_err_v,
		't': obs_err_t,
		'q': obs_err_q
	}).set_index('p').sort_index()
	obs_err = pd.DataFrame(index=df['p'].unique()).sort_index()
	for key in tmp:
		obs_err[key] = np.interp(obs_err.index, tmp.index, tmp[key])
	for i in range(len(df)):
		p = df.loc[i,'p']
		df.loc[i,'u'] += np.random.normal(loc=0, scale=obs_err.loc[p,'u'], size=1)[0]
		df.loc[i,'v'] += np.random.normal(loc=0, scale=obs_err.loc[p,'v'], size=1)[0]
		df.loc[i,'t'] += np.random.normal(loc=0, scale=obs_err.loc[p,'t'], size=1)[0]
		sh0 = df.loc[i, 'sh']
		while True:
			sh = sh0 + np.random.normal(loc=0, scale=obs_err.loc[p,'q'], size=1)[0]
			if sh < 0: continue
			rh = calc.relative_humidity_from_specific_humidity(p * units.Pa, df.loc[i,'t'] * units.K, sh)
			if rh < 0 or rh > 1: continue
			break
		df.loc[i,'sh'] = sh
		df.loc[i,'rh'] = rh.magnitude * 100
	df['d' ] = [x.magnitude for x in calc.dewpoint_from_relative_humidity(list(df['t']) * units.K, list(df['rh']) * units.percent).to(units.K)]
	df['ws'] = calc.wind_speed(list(df['u']) * units.mps, list(df['v']) * units.mps)
	df['wd'] = calc.wind_direction(list(df['u']) * units.mps, list(df['v']) * units.mps)
	return df

def write_sound_littler(file_path, sfc_df, plev_df):
	f = open(file_path, 'w')
	itime = 0
	for time in sfc_df['time'].unique():
		for sid in sfc_df.loc[sfc_df['time']==time]['sid']:
			row = sfc_df.loc[sfc_df['time']==time].loc[sfc_df['sid']==sid].iloc[0]
			plev = list(plev_df.loc[plev_df['time']==time].loc[plev_df['sid']==sid].iterrows())
			# Header
			f.write('{:20.5f}' .format(row['lat']))                           # latitude
			f.write('{:20.5f}' .format(row['lon']))                           # longitude
			f.write('{:40}'    .format(row['sid']))                           # id
			f.write('{:40}'    .format(row['sid']))                           # name
			f.write('{:40}'    .format('FM-35 TEMP'))                         # platform
			f.write('{:40}'    .format('OSSE'))                               # source
			f.write('{:>20.5f}'.format(row['zs']))                            # elevation
			f.write('{:>10d}'  .format(len(plev)))                            # num_vld_fld
			f.write('{:>10d}'  .format(0))                                    # num_error
			f.write('{:>10d}'  .format(0))                                    # num_warning
			f.write('{:>10d}'  .format(itime))                                # seq_num
			f.write('{:>10d}'  .format(0))                                    # num_dups
			f.write('{:>10}'   .format('F'))                                  # is_sound
			f.write('{:>10}'   .format('F'))                                  # bogus
			f.write('{:>10}'   .format('F'))                                  # discard
			f.write('{:>10d}'  .format(-888888))                              # obs_time
			f.write('{:>10d}'  .format(-888888))                              # julian_day
			f.write('{:>20}'   .format(row['time'].strftime('%Y%m%d%H%M%S'))) # date time
			f.write('{:>13.5f}'.format(row['slp']))                           # slp
			f.write('{:>7d}'   .format(0))                                    # slp qc
			f.write('{:>13.5f}'.format(-888888.0))                            # ref_pres
			f.write('{:>7d}'   .format(0))                                    # ref_pres qc
			f.write('{:>13.5f}'.format(-888888.0))                            # ground_t
			f.write('{:>7d}'   .format(0))                                    # ground_t qc
			f.write('{:>13.5f}'.format(-888888.0))                            # sst
			f.write('{:>7d}'   .format(0))                                    # sst qc
			f.write('{:>13.5f}'.format(row['ps']))                            # psfc
			f.write('{:>7d}'   .format(0))                                    # psfc qc
			f.write('{:>13.5f}'.format(-888888.0))                            # precip
			f.write('{:>7d}'   .format(0))                                    # precip qc
			f.write('{:>13.5f}'.format(-888888.0))                            # t_max
			f.write('{:>7d}'   .format(0))                                    # t_max qc
			f.write('{:>13.5f}'.format(-888888.0))                            # t_min
			f.write('{:>7d}'   .format(0))                                    # t_min qc
			f.write('{:>13.5f}'.format(-888888.0))                            # t_max night
			f.write('{:>7d}'   .format(0))                                    # t_max night qc
			f.write('{:>13.5f}'.format(-888888.0))                            # p_tend03
			f.write('{:>7d}'   .format(0))                                    # p_tend03 qc
			f.write('{:>13.5f}'.format(-888888.0))                            # p_tend24
			f.write('{:>7d}'   .format(0))                                    # p_tend24 qc
			f.write('{:>13.5f}'.format(-888888.0))                            # cloud_cover
			f.write('{:>7d}'   .format(0))                                    # cloud_cover qc
			f.write('{:>13.5f}'.format(-888888.0))                            # ceiling
			f.write('{:>7d}'   .format(0))                                    # ceiling qc
			f.write('\n')
			ilev = 0
			for idx, row in plev:
				# Record
				f.write('{:13.5f}'.format(row['p']))                           # pressure (Pa)
				f.write('{:7d}'   .format(0))                                  # pressure qc
				f.write('{:13.5f}'.format(row['h']))                           # height
				f.write('{:7d}'   .format(0))                                  # height qc
				f.write('{:13.5f}'.format(row['t']))                           # temperature (K)
				f.write('{:7d}'   .format(0))                                  # temperature qc
				f.write('{:13.5f}'.format(row['d']))                           # dewpoint (K)
				f.write('{:7d}'   .format(0))                                  # dewpoint qc
				f.write('{:13.5f}'.format(row['ws']))                          # wind speed (m s-1)
				f.write('{:7d}'   .format(0))                                  # wind speed qc
				f.write('{:13.5f}'.format(row['wd']))                          # wind direction (deg)
				f.write('{:7d}'   .format(0))                                  # wind direction qc
				f.write('{:13.5f}'.format(row['u']))                           # wind u component (m s-1)
				f.write('{:7d}'   .format(0))                                  # wind u component qc
				f.write('{:13.5f}'.format(row['v']))                           # wind v component (m s-1)
				f.write('{:7d}'   .format(0))                                  # wind v component qc
				f.write('{:13.5f}'.format(row['rh']))                          # relative humidity (%)
				f.write('{:7d}'   .format(0))                                  # relative humidity qc
				f.write('{:13.5f}'.format(-888888.0))                          # thickness (m)
				f.write('{:7d}'   .format(0))                                  # thickness qc
				f.write('\n')
				ilev += 1
			# End
			f.write('{:13.5f}'.format(-777777.0))
			f.write('{:7d}'   .format(0))
			f.write('{:13.5f}'.format(-777777.0))
			f.write('{:7d}'   .format(0))
			f.write('{:13.5f}'.format(-888888.0))
			f.write('{:7d}'   .format(0))
			f.write('{:13.5f}'.format(-888888.0))
			f.write('{:7d}'   .format(0))
			f.write('{:13.5f}'.format(-888888.0))
			f.write('{:7d}'   .format(0))
			f.write('{:13.5f}'.format(-888888.0))
			f.write('{:7d}'   .format(0))
			f.write('{:13.5f}'.format(-888888.0))
			f.write('{:7d}'   .format(0))
			f.write('{:13.5f}'.format(-888888.0))
			f.write('{:7d}'   .format(0))
			f.write('{:13.5f}'.format(-888888.0))
			f.write('{:7d}'   .format(0))
			f.write('{:13.5f}'.format(-888888.0))
			f.write('{:7d}'   .format(0))
			f.write('\n')
			f.write('{:7d}{:7d}{:7d}\n'.format(ilev, 0, 0))
			itime += 1
	f.close()

def sim_sound(nature_run_root, nature_run_type,
	sid, lon, lat, plev, start_time, end_time,
	sfc_obs_err_u, sfc_obs_err_v, sfc_obs_err_t, sfc_obs_err_q, sfc_obs_err_p,
	plev_obs_err, plev_obs_err_u, plev_obs_err_v, plev_obs_err_t, plev_obs_err_q, ntasks,
	output_littler):
	# Check directory existence.
	if not os.path.isdir(os.path.dirname(os.path.abspath(output_littler))):
		cli.notice(f'Create directory {os.path.dirname(output_littler)}.')
		os.makedirs(os.path.dirname(os.path.abspath(output_littler)))

	# Do jobs in parallel.
	n = len(sid)
	pool = Pool(processes=ntasks)
	if nature_run_type == 'wrf':
		sim_sound_call = sim_sound_wrf
	elif nature_run_type == 'ecwf':
		sim_sound_call = sim_sound_ecwf
	res = pool.starmap(sim_sound_call,
		zip_longest(
			repeat(nature_run_root, n),
			sid, lon, lat,
			repeat(plev, n),
			repeat(start_time, n),
			repeat(end_time, n)
		)
	)

	# Concatenate results.
	sfc_df = []
	plev_df = []
	for i in range(len(res)):
		sfc_df.append(res[i][0])
		plev_df.append(res[i][1])
	sfc_df = pd.concat(sfc_df, ignore_index=True)
	plev_df = pd.concat(plev_df, ignore_index=True)

	# Add observation errors.
	sfc_df = add_sfc_obs_err(sfc_df, sfc_obs_err_u, sfc_obs_err_v, sfc_obs_err_t, sfc_obs_err_q, sfc_obs_err_p)
	plev_df = add_plev_obs_err(plev_df, plev_obs_err, plev_obs_err_u, plev_obs_err_v, plev_obs_err_t, plev_obs_err_q)

	# Write LITTLE_R file.
	write_sound_littler(output_littler, sfc_df, plev_df)

if __name__ == '__main__':
	parser = argparse.ArgumentParser('Simulate sounding observation.')
	parser.add_argument('-i', '--nature-run-root')
	parser.add_argument('-t', '--nature-run-type', choices=('ecwf', 'wrf'), required=True)
	parser.add_argument('-o', '--output-littler', required=True)
	parser.add_argument('--sid', nargs='+')
	parser.add_argument('--lon', nargs='+', type=float)
	parser.add_argument('--lat', nargs='+', type=float)
	parser.add_argument('--plev', help='Pressure levels in hPa', nargs='+', type=float)
	parser.add_argument('--start-time', required=True, type=parse_time)
	parser.add_argument('--end-time', required=True, type=parse_time)
	parser.add_argument('--sfc-obs-err-u', type=float)
	parser.add_argument('--sfc-obs-err-v', type=float)
	parser.add_argument('--sfc-obs-err-t', type=float)
	parser.add_argument('--sfc-obs-err-q', type=float)
	parser.add_argument('--sfc-obs-err-p', type=float)
	parser.add_argument('--plev-obs-err', help='Pressure levels for observation errors', type=float, nargs='+')
	parser.add_argument('--plev-obs-err-u', type=float, nargs='+')
	parser.add_argument('--plev-obs-err-v', type=float, nargs='+')
	parser.add_argument('--plev-obs-err-t', type=float, nargs='+')
	parser.add_argument('--plev-obs-err-q', type=float, nargs='+')
	parser.add_argument('-j', '--json', help='JSON configuration file')
	parser.add_argument('-n', '--ntasks', help='Process number', type=int, default=4)
	args = parser.parse_args()

	if args.json:
		config = parse_config(args.json)
		if not args.nature_run_root:
			args.nature_run_root = config['custom']['nature_run']['dir_pattern']
		if not args.nature_run_type:
			args.nature_run_type = config['custom']['nature_run']['type']
		if not has_key(config, ('custom', 'sim_obs', 'sound', 'obs_err', 'plev')):
			cli.error('There is no pressure level of observation errors in the configure file!')
		loc = np.array(config['custom']['sim_obs']['sound']['location'])
		args.sid = [str(int(x)).zfill(5) for x in np.linspace(1, len(loc), len(loc))]
		args.lon = loc[:,0]
		args.lat = loc[:,1]
		args.sfc_obs_err_u  = config['custom']['sim_obs']['sound']['obs_err']['sfc_u']
		args.sfc_obs_err_v  = config['custom']['sim_obs']['sound']['obs_err']['sfc_v']
		args.sfc_obs_err_t  = config['custom']['sim_obs']['sound']['obs_err']['sfc_t']
		args.sfc_obs_err_q  = config['custom']['sim_obs']['sound']['obs_err']['sfc_q']
		args.sfc_obs_err_p  = config['custom']['sim_obs']['sound']['obs_err']['sfc_p']
		args.plev_obs_err   = config['custom']['sim_obs']['sound']['obs_err']['plev']
		args.plev_obs_err_u = config['custom']['sim_obs']['sound']['obs_err']['plev_u']
		args.plev_obs_err_v = config['custom']['sim_obs']['sound']['obs_err']['plev_v']
		args.plev_obs_err_t = config['custom']['sim_obs']['sound']['obs_err']['plev_t']
		args.plev_obs_err_q = config['custom']['sim_obs']['sound']['obs_err']['plev_q']
	if args.plev_obs_err == None:
		cli.error('Pressure levels of observation errors should be provided!')
	# Convert observation error pressure from hPa to Pa.
	args.plev_obs_err = [p * 100 for p in args.plev_obs_err]

	if args.plev == None:
		args.plev = man_plev
	else:
		args.plev = np.array(args.plev) * 100 * units.Pa

	sim_sound(args.nature_run_root, args.nature_run_type,
		args.sid, args.lon, args.lat, args.plev,
		args.start_time, args.end_time,
		args.sfc_obs_err_u, args.sfc_obs_err_v, args.sfc_obs_err_t, args.sfc_obs_err_q, args.sfc_obs_err_p,
		args.plev_obs_err, args.plev_obs_err_u, args.plev_obs_err_v, args.plev_obs_err_t, args.plev_obs_err_q, args.ntasks,
		args.output_littler)
